Hopping conduction in strong electric fields: 
Negative differential conductivity 
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Effects of strong efectric fiefds on hopping conductivity are studied theoreticaliy. Monte-Cario 
computer simulations show that the analytical theory of Nguyen and Shklovskii [Solid State Com- 
mun. 38, 99 (1981)] provides an accurate description of hopping transport in the limit of very 
high electric fields and low concentrations of charge carriers as compared to the concentration of 
localization sites and also at the relative concentration of carriers equal to 0.5. At intermediate 
concentrations of carriers between 0.1 and 0.5 computer simulations evidence essential deviations 
from the results of the existing analytical theories. 

The theory of Nguyen and Shklovskii also predicts a negative differential hopping conductivity 
at high electric fields. Our numerical calculations confirm this prediction qualitatively. However 
the field dependence of the drift velocity of charge carriers obtained numerically differs essentially 
from the one predicted so far. Analytical theory is further developed so that its agreement with 
numerical results is essentially improved. 

PACS numbers: 72.20.Ht, 72.20.Ee, 72.80.Ng, 72.80.Lc 
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I. INTRODUCTION 



Hopping conduction in solids governed by strong elec- 
tric fields is in the focus of intensive theoretical and ex- 
perimental study since several decades [see, for instance, 
chapter 7 in Ref. [I] and references therein]. In recent 
years, particular interest to this research area has been 
caused by growing device applications of amorphous or- 
ganic and inorganic materials in which the incoherent 
hopping transitions of charge carriers between spatially 
and energetically distributed localized states dominate 
the optoelectronic phenomena [see, for instance, Ref. ,2 
and references therein]. One of the mostly discussed 
topics is whether the differential negative conductivity 
(NDC), i.e., the decreasing conductivity with increasing 
electric field, is possible in the hopping regime. The dis- 
cussion was, to much extent, provoked by the reports 
on the apparent decrease of the drift mobility with rising 
electric field at relatively high temperatures and low field 
strengths in disordered organic materials i^^' 7 -^^ 
This apparent decrease of the mobility with increasing 
electric field was reported to be succeeded by the increase 
of the mobility at higher field strengths. However, the 
self-consistent effective-medium theory for drift and dif- 
fusion at low electric fields^ does not show any decrease 
of the mobility with increasing field. Furthermore, it has 
been shown experimentally 12 and theoretically 13 that the 
apparent decrease of the mobility with rising field at low 
field strengths is an artifact. The experimental data were 
obtained by the time-of-flight technique, in which charge 
carriers are created close to one surface of a sample with a 
given thickness L and the transient time Tt r is measured, 
which is needed for charge carriers to reach the opposite 



surface of the sample at a particular strength of the ap- 
plied electric field F. Then the drift mobility is calculated 
as fi = L/(r tT F). However, at high temperatures and low 
electric fields the current transients in the time-of-flight 
experiments are determined mostly by diffusion of charge 
carriers rather than by their drift. Therefore, using the 
drift formula one strongly overestimates the mobility. It 
is the presence of the field strength in the denominator 
that leads to the apparent "increase" of the mobility at 
decreasing f \ 12 i 13 If one uses at low fields and high tem- 
peratures the diffusion formulas instead of the drift ones, 
then no decrease of the mobility with increasing field can 
be claimed at low electric fields, 12 ! 13 

This result does not exclude, however, the possibil- 
ity of the NDC in the hopping regime. Bottger and 
Bryksin™ and Shklovskii et a L 15 ' 16 have suggested an- 
alytical theories for the mobility and conductivity de- 
creasing with increasing electric field in various disor- 
dered materials. Remarkably, this effect of the negative 
differential conductivity is to be expected at high field 
strengths. This regime succeeds the very strong increase 
of the mobility with rising field ) 15 ! 16 and does not pre- 
cede it at lower fields as claimed on the basis of the drift 
equations, 3 i 4 i 5 i 6 i 7 i 8 i 9 i 10 

The decreasing conductivity with increasing electric 
field at high field strengths has been observed exper- 
imentally for hopping transport in lightly doped and 
weakly compensated crystalline silicon i 17 ! 18 ' 19 The hop- 
ping transport mode in such systems at low electric fields 
had been described theoretically in all detail, 20 which 
made these systems particularly attractive for studying 
the new non-Ohmic effects. Shklovskii et ah 15 i 17 ' 18 i 19 de- 
veloped an analytical theory, which predicted the NDC 
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effect in the lightly doped weakly compensated semi- 
conductors. The experimental observations in lightly 
doped and weakly compensated crystalline silicon appear 
in qualitative agreement with his theoretical predictions. 
Furthermore computer simulations of Levin et ali^i con- 
firmed qualitatively the existence of the NDC effect, 
though no quantitative comparison with the analytical 
theory-- has been attempted. Recent interest in the NDC 
effect has been caused by its importance for construction 
of memory devices. These devices typically contain con- 
ducting particles embedded into a nonconductive mate- 
rial. For such devices, made from both inorgani o 22 ' 23 and 
organi o 24 i 25 ' 26 i 27 ' 28 materials, NDC and switching phe- 
nomena have been reported. Since electrical conduction 
in the materials, which are currently being tried for de- 
vice applications, is dominated by hopping of charge car- 
riers, it is necessary to study the possibility of the NDC 
in this regime in more detail. 

In the present paper we report on the theoretical study 
of hopping transport in high electric fields. In Section 
ILT1 we describe the theoretical model and briefly out- 
line the analytical approach suggested by Nguyen and 
Shklovskii — for the limit of extremely high electric fields. 
In Section[nT]we present our results obtained by straight- 
forward Monte Carlo computer simulations and show the 
range of applicability for the analytical theory of Nguyen 
and Shklovskii. In Sect ion HVl we further develop the ana- 
lytical theory whereby we improve its agreement with the 
results of computer simulations. In particular, the ana- 
lytical dependence of the drift mobility on the concentra- 
tion of charge carriers comes in better agreement with the 
simulation results. Section [V] is dedicated to the NDC. A 
new numerical algorithm has been developed to study the 
NDC effect theoretically. Numerical results obtained in 
the framework of this algorithm confirm qualitatively the 
conclusion of Nguyen and Shklovskii on the possibility of 
the NDC in the hopping regime. However, the field de- 
pendence of the drift velocity of charge carriers obtained 
numerically differs essentially from the one predicted so 
far^ We suggest in Section [V] a further development of 
the analytical theory, improving essentially its agreement 
with numerical results. Concluding remarks are gathered 
in Section IVT1 



II. MODEL AND THEORETICAL 
BACKGROUND 

Aiming to clarify whether the NDC effect is inherent 
for the hopping transport regime, we consider first, fol- 
lowing Nguyen and Shklovskii, 15 the simplest possible 
model — a three-dimensional array of isoenergetic sites 
with a random spatial distribution with the concentra- 
tion TV. Each site can be either empty or occupied by a 
single electron. Energies of electrons are equal on all sites 
so that no energy disorder and no electron-electron inter- 
actions between different sites are taken into account. 
Only in the final part of Section [V] we study the ef- 





FIG. 1: The shape of optimal traps for the infinite (a) and a 
finite (b) electric fields. The dotted path in (a) is forbidden 
at infinite fields but provides an escape route at finite fields 



feet of the energy disorder on the NDC. An electric field 
F = (— F, 0, 0) is put along the negative direction of the 
axis X, so that the drift velocity of the negatively charged 
electrons is directed along the X axis. Conduction takes 
place due to tunnelling hops of electrons between the lo- 
calization sites. The rate r^- for an electron hop from 
site i to site j is determined as 



Fij = T exp 



a 



eF{ Xj 



kT 



rii(l-nj), (1) 



where is the distance between the sites, a is the lo- 
calization length, e is the elementary charge, k is the 
Boltzmann constant, T is the temperature, and n^rij 



are the occupation factors of the sites, 



ni,rij 



e{0;l». 



The function / is related to energy the gain or the energy 
loss during the jump: 



/(«) 



I, ifa>0, 
exp(a), if a < 0. 



infinite electric field, the factor 
reduces to the Heaviside's function 



In the limit of 

/ [eF( Xj - Xi)/kT 
9{xj - Xi). 

Below we will assume that To = e = 1. As a measure 
of length, we introduce the typical distance between the 
neighboring sites R = TV" 1 / 3 . 

Using this simple model, Nguyen and Shklovskii 15 have 
shown analytically that the effect of the NDC is inherent 
for hopping transport. Let us consider briefly their ar- 
guments starting from the case of infinitely high fields F 
and extremely small electron concentrations n e . Under 
such circumstances each electron can be treated inde- 
pendently from the others and electrons can move only 
toward the increasing values of their x coordinate. In 
Fig. [T] this is the direction to the right. At each jump, 
an electron moves along the axis A to a distance ~ R, 
so that its drift velocity can be estimated as v ~ R/f, 
where r is an average time between jumps (dwell time). 
A dwell time Tj for hopping from the site i is of the or- 
der of exp(2r^/a), where r, is the distance from site i to 
its nearest neighbor "to the right", i.e. with co-ordinate 
x larger than Xi. In other words, r» is the maximum 
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radius of a hemisphere centered at the site i that does 
not contain any other sites (Fig. DJa)). If is much 
larger than R (the typical distance between the neighbor- 
ing sites), such an empty hemisphere can be considered 
as a trap for electrons. The contribution of traps with 
radii in the range [r, r + dr] to the average dwell time 
t is proportional to r(r) = exp(2r/a), and also to the 
probability of the corresponding configuration of sites, 
p(r)dr = 2nNr 2 exp(-2?riVr 3 /3)dr: 



r(r)p(r)dr 



2nNr z exp 

a 



2r 2ttN 



dr. 



This integral is easy to evaluate, taking into account 
that the integrand has a sharp maximum at r = r m = 
l/y/nNa. Consequently one obtains for the current den- 
sity j = n e v 



n e R 



JF- 



T 



n e (a 3 R) 1/4 exp 



30? 



R 



Therefore, one can conclude that in the limit F — > 
oo,n e — > the current is determined by hemispherical 
traps (Fig.[T^) with an "optimal" radius r m = 1/yTvNa. 

In the case of finite electric fields, a hemispherical trap 
is not an efficient one, because an electron has a possi- 
bility to move in the energetically unfavored directions, 
and thus to escape the trap (for example, along the 
dotted arrows in Fig. [Ua)). According to Nguyen and 
Shklovskii, 15 an "optimal" trap for an electron in large 
though finite electric fields F consists of a hemisphere to 
the right and of a cone to the left of the site on which 
an electron is captured, with a chain of sites along the X 
axis that provides an easy path into the trap (Fig. [TJd). 
The height h of the cone is chosen so that it is equally 
hard to escape the trap in all directions taking the chain 
along the X axes into account: h — 2rkT / Fa. Therefore, 
the smaller is the field, the larger is the volume of a trap 
with the same dwell time, and consequently the smaller is 
the probability p(r) of finding such a trap. It means that 
the average dwell time r = J Q r(r)p(r)dr decreases with 
decreasing field strength, and concomitantly the current 
density j ~ n e i?/r increases with decreasing field. This 
is the essence of the physical mechanism that causes the 
NDC effect.— To obtain an expression for the current 
density, one can substitute the volume of the trap shown 
in Fig. [lb, Vtrap = (1 + kT/Fa)2nr 3 /3, instead of the 
hemispherical trap volume, 27rr 3 /3, into the integral |2]). 
The result reads 



n e (a 3 R)^ 



exp 



30F 



R 



kT 
~Fa 



(4) 

This is the mathematical expression for the NDC. The 
approach leading to this expression is applicable only for 
fields F 3> kT/R. In smaller fields, the assumption that 
almost every jump is directed along the axis X is vio- 



lated. Therefore one should expect that Eq. ([4]) overes- 
timates the current density for F ~ kT/R. 

Equations ((3]) and |(3} are valid only if the concentra- 
tion of electrons n e is small as compared to the concentra- 
tion of "optimal" traps, n m = N exp[— NVt Tap (r m )]. In 
the opposite case, n m <C n e <C N, the "optimal" traps 
are almost always occupied and play a negligible role. In 
such a case the most important traps, which determine 
the drift velocity of electrons, are the ones whose con- 
centration is equal to n e . One can estimate the electron 
drift velocity as v — l/T(r n )n e S , where r n is a radius of 
the most important traps, r(r n ) is their dwell time, S is 
their capture crosssection. Assuming S ~ R 2 = iV~ 2 / 3 , 
one obtains for the current flow j — n e v ~ iV 2 / 3 r(r I1 ) _1 . 
For an infinitely large field, the radius r n is defined via 



7Vexp(-W trap (r„)) = TV exp 



/ ., N \ 1/3 

that gives r n = R I ^ log j . Consequently, the cur- 
rent density isi^ 



2ttA 



AT 2 / 3 

T ( r n) 



Nz exp 



2 / 3 , N 

- h^-log — 

a \ 2ir rip 



The corresponding expression for the concentration 
range n m -C n e <€. N in the case of finite electric fields 
was also obtained by Nguyen and Shklovskii [see Eq. (11) 
inRef.GJ]. 

The case of almost filled sites, n e i=s N, is similar to the 
case of almost empty sites, n e « due to electron-hole 
symmetry. The current density is a symmetrical function 
of the electron concentration: j(n e ) — j(N — n e ). 

Nguyen and Shklovskii^ also emphasized that a spe- 
cial consideration is needed for the case of half-filled sys- 
tem, n e = N/2. They have shown that the concept of 
directed percolation can be used to obtain the current 
density at infinitely high electric fields. In the half-filled 
system the trapping of electrons does not play any role, 
because (due to the electron-hole symmetry) it does not 
change the electron concentration on the infinite cluster 
which is responsible for the current. Current is deter- 
mined by electron jumps to distances d E [r^rf + a/2], 
where r c \ is the percolation threshold of a directed per- 
colation problem. The number of pairs of sites with dis- 
tances d S [rf,rf + a/2] in the infinite cluster per unit 
area is \/L\, where L± — R(2rf/a) lJ± is a transversal 
correlation length of the percolation cluster, and is± is a 
critical indexJ^ The current density is equal to^ 



1 



n„=l/2 



Llr{ri) 



= 7V3 



a 

~2~A 



2l/J_ 



exp 



2r° 



(6) 

Nguyen and Shklovskii^ have also obtained the value of 
the percolation threshold r d c = (0.93 ± 0.01) R and that 
of the correlation length index v = 1.2 ± 0.1. 
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The above arguments of Nguyen and Shklovskiii^ pro- 
vide an analytical theory of non-Ohmic hopping conduc- 
tion, based on the concept of the trapping-determined 
transport. The theory is valid for the case of large elec- 
tric fields in two concentration ranges: n e -C n m and 
n m < n e < JV. Most remarkably, this theory predicts 
the effect of the NDC. Also a theory for the case of the 
half- filled system (n e — N/2) for infinitely high electric 
fields (F — > oo) has been suggested based on the directed- 
percolation-approachJ^ 

Below we present our numerical study of the field- 
dependent hopping conductivity. It shows the range 
of validity for the analytical theory of Nguyen and 
Shklovskii. 15 Furthermore, the analytical theory is devel- 
oped below in order to improve the agreement between 
the analytical and numerical results. 



III. MONTE CARLO SIMULATIONS FOR 
INFINITELY HIGH FIELDS 

In order to calculate the electron drift velocity and the 
current density at high fields, we used a Monte Carlo 
approach. In the limit of infinitely high fields the di- 
rection of the electron motion is prescribed. Therefore 
it was possible to simulate by a Monte Carlo algorithm 
the motion of an electron in an infinite medium along 
the field direction and therefore to avoid any size effects. 
Without loosing generality one can restrict the maximal 
length of electron transitions involved into the algorithm 
by a reasonably large value d max . In order to simulate 
the fc-th Monte Carlo step in the electron motion, one 
has to store information only about sites inside a layer 
Xk < x < Xk + dmax, where Xk is the electron coordi- 
nate before the A:-th step. We have chosen d max = 3R, 
which provides a possibility to hop to 27rd^ lax A r /3 ~ 57 
neighbors in average. For all sets of parameters used in 
the simulation the size of the optimal trap r m considered 
by Nguyen and Shklovskii was essentially less than d max . 
Therefore, the restriction imposed by d max did not lead 
to any loss of generality. Before making the next step, 
the computer can forget all the information about sites 
in the layer Xk < x < Xk+i, but it has to get information 
about new sites in a layer Xk + d max < x < Xk+i + d max . 
As these "new" sites did not affect the calculation at all 
previous steps, they can be created at random. There- 
fore, each Monte Carlo step includes not only the choice 
of a jump, but also a generation of some "new" sites 
and deleting some "old" sites. To make their number 
finite, one should restrict the system size in the direc- 
tions perpendicular to the field, i.e. in the plane YZ. A 
calculation domain < y < 120R, < z < 120R with pe- 
riodical boundary conditions in the plane YZ was used. 
The motion of a single electron was simulated within the 
described algorithm in order to evaluate the drift velocity 
in the limit n e — > 0. 

For finite electron concentrations, we perform simu- 
lations in a cubic domain with size 60R x 6QR x 60i? 



(a) 








If 






Eq. (3) 

Eq. (11) . 
simulation • 



0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
a(R) 



co 0.4 

-z. 

> 

> 0.3 





(b) 




**** • 










< 






Eq. (11) 
simulation i — • — i 
i i i i i i i i 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
a(R) 

FIG. 2: (a) Drift velocity v in the limit of infinitely large 
electric field and small electron concentration, as a function 
of localization length a; (b) the same data divided by the 
velocity value predicted in Ref. Il5l (Eq. ©). 



and with periodic boundary conditions for all three axes. 
The rates of all possible jumps are calculated before 
starting the Monte Carlo steps, but without the factor 
rii(l — rij) related to occupation. This factor determines 
which jumps are allowed and which are forbidden. Infor- 
mation about allowed and forbidden jumps is updated at 
each step. We used a binary-tree data structure for stor- 
ing the jump rates that gives the possibility to "switch 
on" and "off" jumps efficiently. 

A routine Monte Carlo procedure has been used. In 
each Monte Carlo step the final site for electron hops 
was calculated via the probabilities proportional to the 
hopping rates to different sites and the time At spent to 
jump was calculated via the reciprocal of the sum of rates 
of all possible jumps. Hops from an occupied site were 
possible to any empty one in the direction of increasing 
coordinate x with the restriction that the hop distance 
is less than <i max = 3i?. At each hopping event the in- 
crement Ax in electron x-coordinate is calculated. An 
outcome of the simulation is either an average velocity of 
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an electron, 

in the case of single electron hopping, or a flow of elec- 
trons, 

in the case of finite electron concentrations (where O is 
the volume of the calculation domain). The summation 
was carried out over all sequential Monte Carlo steps. 
For simulations of the behavior of a single electron in an 
empty system, 10 7 Monte Carlo steps were used for a ^ 
0.2R, 10 s steps for a = O.lOi? and 0.15R and 10 9 steps 
for a = 0.07-R. As a result, for a 0.157? convergence 
was not worse than 1%. 

For finite electron concentration, 5 • 10 7 Monte Carlo 
steps were used, this gave a convergence not worse than 
1% for a given realization. At a ^ OAR, there were 
sometimes essential differences between current densities 
in different realizations. The scatter is shown by error 
bars in the figures. 

Simulation results for the electron drift velocity v = 
j /n e in the limit n e — > are shown in Fig. [U(a) by dots 
as a function of the localization length. The analyti- 
cal result of Nguyen and Shklovskii (Eq. ©) is shown 
by the solid line. One can see that Eq. ([3]) correctly 
describes the dependence of the drift velocity on the lo- 
calization length and, furthermore, it correctly estimates 
the magnitude of the velocity. The concept of Nguyen 
and Shklovskii on the hopping drift velocity controlled 
by hemispherical traps is herewith confirmed. However, 
there is some deviation of the simulation data from the 
analytical results. To make this deviation more trans- 
parent, we plot the ratio of the simulated drift velocity 
to its analytical prediction (Eq. ©) in Fig. [D(b). It is 
seen that Eq. ([3]) overestimates the electron velocity by 
a factor of two to five. In Section HvK . some reasons for 
this mismatch will be considered. The analytical theory 
is further developed there to give a better agreement with 
the simulation data. The result of the improved theory 
for the drift velocity (Eq. (fTTj) ') is also shown in Fig. [2]by 
the dashed line. 

The dependence of the current density on the electron 
concentration is shown in Figs. [3] and [U Fig.[3]shows this 
dependence in a wide concentration range, in comparison 
with the analytical results for small (n e <C n m , Eq. 
solid line) and intermediate (n m <C n e -C N, Eq. |5]), 
dashed line) concentrations. One can see that the sim- 
ulated concentration dependence can be roughly divided 
into three parts: for very low concentrations (n e < n m ) 
the dependence is linear^ in accordance with Eq. ([3]); 
then, for n m < n e < 0.03N, it becomes superlinear, as 
described by Eq. ([5]); and finally, for n e > 0.03A, this 
dependence is sublinear and is not described by the the- 
ory based on the transport controlled by traps. In Sec- 
tion [IV}3, we will present an analytical approach valid for 
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FIG. 3: Current density j as a function of electron concentra- 
tion n e for values of the localization length 0.2R, 0.1R, and 
0.07i? (from top to bottom). The electric field is infinitely 
large. 



the range of parameters covering the ranges of applica- 
bility of Eqs. © and ([5]). The result of this developed 
approach is Eq. (fT2)) shown by dashed lines in Fig. [3^ 
One can see that it provides an accurate description of 
the current density for any concentration less than 0.03iV. 

For n e > 0.03^, the simulated values of the current 
density are smaller than those predicted by the analyt- 
ical theory due to the following reason. At sufficiently 
large electron concentrations, the conducting paths are 
not almost empty, as is assumed in the theory. More- 
over, there are "bottlenecks" for the current, where the 
electron concentration is much larger than the mean con- 
centration n e . In these places, the factor of (1 — rij) in 
Eq. (fTJ) turns out to be important, and due to this factor 
the current density is suppressed. 

In Fig.[H the simulation results are shown for the whole 
range of carrier concentrations. For convenience, all val- 
ues of current density are divided by the maximum value 
for the given localization length. For large localization 
lengths (a > 0.15i?), the concentration dependence of 
the current density j obeys approximately a parabolic 
law: j(n e ) ~ n e (N — n e ) at concentrations in the vicin- 
ity of the half filling. One can interpret this behavior in 
terms of the hopping rates, namely, by substituting the 
mean occupancy n e /N instead of rii and rij into Eq. (fT]). 
Concomitantly, one obtains that the contribution of each 
pair of sites is proportional to n e (N — n e ). The same 
concentration dependence is expected then for the cur- 
rent density. 

Fig. [5] shows the simulated dependence of the current 
density on the localization length (dots) in comparison 
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FIG. 4: Current density j as a function of electron concentra- 
tion n e , normalized on its maximum value j(N/2), at different 
localization lengths. The electric field is infinitely large. 




FIG. 5: Current density j as a function of localization length 
a for the half-filled system (n e = N/2). The electric field is 
infinitely large. 



with the analytical theory based on the concept of di- 
rected percolation (Eq. ([5]), solid line) for n e = N/2. 
Apparently, the theory of Nguyen and Shklovskii— cor- 
rectly describes this dependence within the range of cur- 
rent densities of almost 15 orders of magnitude. However, 
the theory underestimates the magnitude of the current 
density by approximately a factor of 30. Further research 
is necessary to clarify the reasons of this discrepancy. 



IV. ANALYTICAL THEORY FOR INFINITELY 
HIGH FIELDS 

Our numerical studies show that although the ana- 
lytical description of hopping conduction in very strong 
electric field by Nguyen and Shklovskii is qualitatively 
correct the quantitative predictions differ sometimes by 



more than an order of magnitude from the numerical re- 
sults. In this Section, we show how to improve the accu- 
racy of the analytical theory. 



Limit of n e 







For low electron concentrations n e <C n m , where n m — 
iVexp[— ^-(R/ira) 3 / 2 ] is the concentration of the "opti- 
mal" traps, the prediction of Ref.[l5| for the electron drift 
velocity is expressed by Eq. Now we discuss several 
corrections to this equation. 

1) There is a numeric factor of (47T) 1 / 4 w 1.88 m r, 
arising from the evaluation of the integral ((2|) that should 
be taken into account. It gives a factor of (47r) -1 / 4 for 
the drift velocity. 

2) The mean electron displacement along the X axis, 
{Ax), is taken equal to R in Ref. [l5|. We performed 
Monte Carlo calculations for (Ax) as a function of the 
localization length a and get the following fitting expres- 
sion: 



(Ax) = R (0.385 + 0.45a 2 /i? 2 



(7) 



(the accuracy of fitting is not worse than 0.3 % in the 
range 0.05 < a/R < 0.2). Therefore, the drift veloc- 
ity v = (Ax) /t gets an additional factor equal approxi- 
mately to 0.5. 

3) The dwell time r(r) of a trap with a radius r is 
in fact somewhat less than the value exp(2r/a) used in 
Ref. [HI because an electron can escape the trap by mov- 
ing not only to the nearest site to the right, but also to 
a more distant site. A contribution Ti of these distant 
sites to the escaping rate is 



I\ = J (T 2rila 2KNr\dr 1 = e - 2r/a irNa 



r +ar+ — 



Then, the dwell time r(r) is a reciprocal of the sum Tq + 
Ti, where Tq = exp(— 2r/a) is the rate of a jump to the 
nearest neighbor: 



r(r) 



1 



exp(2r/a) 



r, 



1 + irNa(r 2 + ar + a 2 /2) ' 



(8) 



For r — r m = (TrNa)^ 1 ^ 2 , r(r) is approximately half 
of exp(2r/a), that results in a factor of two in the drift 
velocity. 

4) The geometrical crosssections of larger traps have 
larger capture crosssections for electrons than the smaller 
ones. This results in different probabilities for carriers to 
be captured by traps with different radii. The probability 
p(r)dr that the next visited site will be a trap with radius 
in the range (r, r + dr) is 



p(r)dr 



5(r) 
(S) 



p(r)dr, 
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a=0.05 R 
a=0.07 R 
a=0.10 R 
a=0.15 R 
a=0.20 R 
a=0.25 R 
fitting: 0.81 + 0.36 r 2 /R 2 
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FIG. 6: Relative capture crosssection as a function of squared 
trap radius. 



where p(r) = 2irNr 2 exp(-2?riVr 3 /3), S(r) is a cap- 
ture crosssection of a trap with radius r, and (S) 
| S(r)p(r)dr is a mean crosssection. Below we will use 
a notation SVcifV) for a "relative crosssection" S(r)/(S). 
Then, instead of Eq. @ we get 



oo 

t = j V(r)S r 



ei{r)p(r)dr. 



(9) 



We calculated the relative crosssections with Monte 
Carlo method as ratios -/Vt^r, r+Ar]/(Njp(r)Ar), where 
^Vtr [r, r+Ar] is a number of traps with radii in the speci- 
fied range visited by an electron, and Nj is a total number 
of electron jumps. We used Nj = 10 8 and Ar — O.Oli?. 
The results are presented in Fig. [S] The relative crosssec- 
tion is almost independent of the localization length for 
r > 0.3R. Its dependence on the trap radius is described 
by the quadratic function: 



S le i(r) = 0.81 + 0.36 r 2 /R 2 



(10) 



For the "optimal" traps with r = r m = (nNa) x l 2 we 



get S Te i(r m ) ~ a in the limit 
Eq. © , it results in a factor of ~ 



-> 0. According to 
for the mean dwell 



time t and, consequently, in a factor of ~ a for the drift 
velocity. 

Now we can improve Eq. ([3]) of Nguyen and Shklovskii, 
starting from Eq. ©. Since the integrand has a sharp 
maximum at r rn = (irNa) -1 / 2 , we can estimate the inte- 
gral approximately as 

t « T(r m )S Te i(r m )p(r m ) ( 7 rfl 3 a/4) 1 / 4 iV- 1 . 



Then, using Eqs. Q, ©, and (TO}, we get the following 
expression for the drift velocity v = (Ax) jf: 



(0.85 + 0.45^) [l + nNa ^ + ^+^ 



(4tt) 1 /4 (0.81 + 0.36 



7T( ! / 



Here vns = i/*n. e is the drift velocity corresponding to 
Eq. (|3|) . Finally, the latter expression can be fitted (with 
accuracy of about 3 % for a < 0.2R) by a simple formula, 



a 
R 



(l.4 + 2.1e- 10a / 7? 



(11) 



This expression is to be considered as a corrected analyt- 
ical form for the drift velocity at infinitely high fields in 
the limit of small electron concentration. 

A comparison of Eq. (JTTJ) with the values of the drift 
velocity obtained by the Monte Carlo method is shown 
in Fig. [2 The difference between the analytical and sim- 
ulated results does not exceed 20 %. We believe that the 
main source of this small difference is some inaccuracy 
in determining r(r) by Eq. l[8]). In fact, for a given trap 
radius there is some distribution of the dwell times. The 
quantity r(r) contributing to Eq. @ is the mean dwell 
time for radius r. However, Eq. |(5J) gives the reciprocal 
value of the mean escaping rate that is slightly smaller 
than r(r). For this reason, Eq. (JTTJ) can slightly overes- 
timate the drift velocity. 



B. Finite electron concentration 

Let us now try to improve the analytical approach at 
finite, though small electron concentration, n e <C N/2. 
In this case, electron flow can be considered as a ho- 
mogeneous one on the scale of distances between the 
traps that determine the transport. Hence one can ex- 
press the frequency z/;„ of electron capture by a trap as 
Vin = jS(l — n), where j is the current density, S is 
the trap capture crosssection, and n is its mean occu- 
pancy. Under the steady-state conditions, Vi n = v ou t, 
where v out — h/r is a frequency of emission of electrons 
from the trap, r is a dwell time. From this equation one 
can get n: 

1 



1 + 0'Sr)- 1 



Since in a snapshot of the system almost all electrons are 
captured by traps, the total electron concentration n e is 



n(r)p(r)dr 



p(r)dr 



l + (j5(r)r(r))-i' 



(12) 



The dwell time r(r) can be estimated by Eq. ([5]). In or- 
der to find the crosssection S(r), one should note that 
the mean crosssection (S) is equal to l/N(Ax). Conse- 
quently, 



S(r) 



SVdM 2 0.81 i? 2 + 0.36 r 2 
N(Ax) ~ 0.385 R 2 + 0.45 a 2 ' 



(13) 



Equation (|12|) with r(r) and S(r) determined by 
Eqs. ([5]) and (|13p . respectively, gives a functional depen- 
dence between the electron concentration and the current 
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density for any n e <C N/2. Fig. [3] evidences a good agree- 
ment between Eq. (Ti"2"|) and the Monte Carlo calculations 
for n e < 0.03 N. 

Although there is probably no simple way to resolve 
Eq. (fT2)) with respect to j analytically in the general case, 
it is possible to simplify this equation in some limiting 
cases. For small n e and j (n e -C n m ), the unity term 
in the denominator of Eq. (|12[) can be dropped, and we 
get n e — j/v, where the drift velocity v = (Ax)/? is 
determined by Eq. (fTTjl . In the opposite limit (n e ^> n m ), 
one can evaluate Eq. (fT2|) as 



27rJVr3 \ 



p(r)dr = iVexp 



where a cutting parameter r n is given by condition 
jS(r n )T(r n ) = 1. Therefore, 



1 



J 



S(r n )T(r n ) 



with 



«lflog~ 

Z7T n e 



1/3 



(14) 



(15) 



Equation (| 14[> . with parameters determined by Eqs. (J8j> , 
(TL3"]) and (I15[) is the improved version of Eq. © by 
Nguyen and Shklovskii for the concentration range n m <C 
n e < N/2. 



V. HOPPING TRANSPORT AT FINITE 
ELECTRIC FIELDS 

So far we have considered the limiting case of infinitely 
high electric fields. Let us now turn to the field depen- 
dence of the charge carriers velocity in order to reveal 
the possibility of the negative differential conductivity 
predicted by Nguyen and Shklovskii. 15 Eq. (g]) predicts 
a decreasing drift velocity with increasing electric field, 
provided the field is strong enough. On the other hand, 
for very small fields, Ohmic transport can be expected, 
i.e., the drift velocity should depend linearly on the field. 
In order to simulate hopping transport at finite electric 
fields, we solved a system of balance equations instead of 
using a direct MC simulation. In the following subsec- 
tion A we describe the details of the numerical procedure, 
while the results are presented in subsection B. 



A. Balance equation method 

We consider a cubic system (side length L) with ran- 
domly placed sites. Periodic boundary conditions are 
used in all directions. The balance equation for the occu- 
pation probability pi of a site i has the for m 29 i 30 ' 31 i 32 i 33 i 34 

^PiTijil-pj) = ^pjFjiil-pi). (16) 



If all occupation probabilities pi are small, i.e. the 
charge carrier concentration is low, the balance equation 
can be linearized: 



ji- 



(17) 



These equations are solved by defining 











r 2 i 


r 3 i • 


■■) 
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-r 2 


T32 ' 
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P3 
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Tl3 
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_r 3 . 






\ J 




I '■ 









(18) 



where T{ = J2j=£i ^ij ^ s the rate of jumping out of site 
i. The equation is then Mp = 0, which we solve nu- 
merically. The matrix M defined in this way is singular, 
which makes a direct solution rather difficult. By replac- 
ing one of the balance equations with the normalization 



E 



Pi 



i. 



(19) 



the matrix becomes nonsingular, and the solution can 
be obtained more efficiently. Additionally, the solution 
obtained in this way is correctly normalized. After this 
replacement, the equation has the form: 



111 

ri2 — r 2 
L13 r 2 3 — L3 



V 



\ 


(pi) 








P2 









P3 







J 


V J 




W 



(20) 



As in Section [TO] the rates for jumps longer than d max 
are assumed to be zero. Hence, it is efficient to use a 
sparse storage scheme for the matrix, where only the non- 
zero elements are stored. We obtained the most accurate 
results in the shortest time by solving the equation by LU 
factorization (in Matlab or Octave with the \ operator). 
This method demands much memory, and did not work 
for L above about 22 R on a 32-bit computer. 

When the steady-state occupation probabilities are 
known, the average velocity of a charge carrier along the 
field direction is given by 

(v x ) = X ~ x 0i ( 21 ) 

and the mobility is then pL = 



B. Field dependence of the current density 

The simulated dependence of the drift velocity v on 
the electric field F is presented in Fig. [7] Simulations are 
performed for 20 3 sites in a cubic domain with the size 
L = 20 R, in the limit of infinitely small electron concen- 
tration. Different symbols refer to different localization 
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FIG. 7: Drift velocity as a function of the electric field, for 
different localization lengths. The system size L is 20 R. 




0.25 



FIG. 8: The field F m corresponding to the maximum of the 
drift velocity as a function of the localization length a. The 
linear fit is given by F m R/kT = YJa/R - 0.2. 



lengths and/or different realizations of the distribution 
of sites in the domain. The size of the simulated system 
was 10 times larger than that in the simulations of Levin 
et al.)2I whose computer simulations for the first time 
confirmed the existence of the NDC effect for hopping 
transport. 

In the limit of small electric field, FR/kT <C 1, sim- 
ulations show an Ohmic conductivity, i.e., v is propor- 
tional to F, in accordance to the Miller-Abrahams con- 
cept of the resistance networki 20 i 35 With increasing field, 
the drift velocity reaches a maximum value. The field 
strength F m corresponding to the maximum of the ve- 
locity appears to be nearly proportional to the localiza- 
tion length a within the range 0.08 R < a < 0.2 R (see 
Fig. [8]). At field strengths F > F m the NDC appears, i. e. 
the drift velocity drops with increasing field. Simulations 
show the presence of the NDC for localization lengths up 
to 0.3 R; when the localization length is decreased, the 
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FIG. 9: Drift velocity as a function of the electric field for 
different localization lengths. The curves show Eq. ((4]) and 
Eq. (|24p . scaled to approach Eq. in the limit of large 

fields. The system size L is 20 R. 



NDC effect becomes more pronounced. 

Fig. [9] shows the comparison between the simulation 
results (symbols) and the predictions of Nguyen and 
Shklovskii 15 (Eq. ([4]), dashed lines). For better agree- 
ment between the theory and the simulation at F — > oo, 
we take into account the i^-independent correction (fTT|) 
to Eq. (H]). Some discrepancies between the simulated 
and predicted drift velocities remain at large fields for 
a = 0.10 R and a — 0.12 R. We believe that these dis- 
crepancies are due to the small size of the simulated sys- 
tem. In fact, the simulated system must be large enough 
to contain a reasonable number of optimal traps. The 
concentration of optimal traps decreases sharply with de- 
creasing localization length, so that at smaller localiza- 
tion lengths larger systems are needed. Thus, for small 
localization length (0.10 R and 0.15 R), only the shape of 
the simulated filed dependence should be taken as repre- 
sentative, but not the values of the calculated velocities 
themselves. 

The range of applicability of Eq. ^ is determined by 
the condition FR/kT ^> 1. One can see nevertheless that 
even within this range the field dependence of the velocity 
is much weaker than the one predicted by Eq. (j4|). This 
result forced us to consider another possible optimal trap 
shape for the case of a finite field as compared to the one 
considered in Ref. O 

The essential feature of the optimal trap proposed by 
Nguyen and Shklovskii (Fig.Q})) is the chain of sites along 
the axis of the cone. This chain was introduced in order 
to provide an easy path for an electron into a trap. The 
chain affects the trap shape and volume, as it serves also 
as a channel for escaping of an electron from the trap. We 
suggest that at moderate localization lengths (a ~ 0.1 R) 
traps without such a chain can also play a significant role. 
Our next aim is to consider the shape of traps without 
a chain of sites and to estimate their influence on the 
electron conduction. 
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FIG. 10: Geometry of the optimal trap at finite electric fields 
without a chain of sites leading into the trap. 



A sketch of such a trap is shown in Fig. [TO] Its shape is 
defined by the condition that the rate of jumping from the 
central site to any point of the trap's surface is the same. 
From Eq. ([T]) one can see that in the positive direction 
along the axis X the trap is bounded by a hemisphere, 
and in the negative direction — by a surface defined by an 
equation 



2r 
a 



(22) 



where r is the radius of the hemisphere, and the origin 
is placed at the central site of the trap. The surface 
determined by Eq. (|22|) is a quadric surface (an ellipsoid, 
a paraboloid, or a hyperboloid, depending on the values 
of parameters). The volume of the trap, 



Krap is 



Krapto 



2?rr 3 



2(c+l) 2 



(23) 



where c = Fa/2kT. 

To obtain the drift velocity v (or the current density 
j = n e v) in the assumption that the most important 
traps are those shown in Fig. [TOl one can proceed in the 
same way as the one applied in Section [ill to get Eq. ([4]); 
the only difference is using V^ rap (r) instead of Nguyen 

and Shklovskii's trap volume Vt rap (r) = (l + 
The result is 



, 3 „ N i 4 ( R 

jn e ->o ^ n e (a R)* exp - — p= — 

\ a 

For the optimal trap radius one gets 



c + 2 
2(c+l) 2 



(24) 



1 



yJnNa 



1 



2(c+l) s 



-1/2 



Since V t ' rap (r) 



< VtrapWj the new trap shape gives a 
weaker field dependence of the drift velocity, and a better 
agreement with the data from simulations, as one can 
see in Fig. |H1 However, the simulated field dependence 
appears even weaker than the one expressed by Eq. ([24]) . 
It leads to the assumption that an actual optimal trap 
has a shape different from both Fig. []J> and Fig. [TQ1 and 
hence has a different volume. 
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FIG. 11: Time-averaged density of sites around the charge 
carrier for different fields at the localization length a — 0.15 R. 
The position of a charge carrier (at the origin) is pointed out 
by a cross. Boundaries of optimal traps predicted in Ref. [l5| 
(see Fig. [1} are depicted by solid lines, the boundary of the 
trap sketched in Fig. [10] by the dotted line. Spatial coordi- 
nates are in units of R. The value 2.36fcT for FR corresponds 
to the maximum of the drift velocity. 



To further investigate the shape of the most efficient 
traps, we collect information about the trap shape from 
the simulations. Fig. QT] show the time-average of the 
density of sites around the charge carrier. To calculate 
this density p(r), the space was divided into small ele- 
ments of equal volume AV; then p(r) was evaluated as a 
sum over pairs of sites: 



p(r) 



1 

AV 



where pi is an occupation probability of the i-th site; 
Xij{ r ) = 1 if the vector (r^ — r,) belongs to the same spa- 
tial element as the vector r; otherwise Xij ( r ) = 0. Finally, 
values of p(r) were averaged over several realizations of 
site distributions. 

Since the carrier spends the most time in the efficient 
traps, the density distribution p(r) directly reflects the 
shape of these traps. At high fields (Fig. [TTk ) the hemi- 
spherical shape and the size of the trap are in an excel- 
lent agreement with the Nguyen and Shklovskii's theory. 
However, at moderate fields, in the region of the NDC 
(Fig. HHb)), neither Fig. C^b) nor Fig. [TO] describe the 
simulated optimal trap. The optimal trap consists in 
such a case of a hemisphere in the spatial region x > 0, 
and of a toroidal "barrier" in the region x < 0, adjoining 
to a periphery of the hemisphere. The volume of the op- 
timal trap turns out to be smaller than the one predicted 



11 



0.0 













o/kT = 


0.0 


1 












o/kT = 


0.2 




f 1 










O/ft 1 — 


n A 












o/kT = 


8 














o/kT = 


i!o 












o/kT = 


1 5 




□ ■ 


* 








o/kT = 


2.0 






□ 
■ 








o/kT = 


3.0 




mm 








o/kT = 


4.0 




n o ° 
Si ° 

x " o 










o/kT = 


5 




G 


s 






o/kT = 


io!o 




% o 




o 












° • 


• 


• 


• 


s 

A 

A 








1°** 


A 


A 


A 

A 






! 

A 

V 




A 


A 

— * 


V 

— 4- — 


V 

— t— 






T 


5 




10 


15 


20 


25 




30 








FR (kT) 









35 



site (with a Gaussian distribution) decreases the drift 
velocity and it also decreases the height of the peak of 
the velocity as a function of the electric field as compared 
to systems with only spatial disorder. This weakening of 
the NDC effect with the increase of the energetic disorder 
(or with the decrease of temperature) is in agreement 
with experimental observations! 18 > 19 Generally, the 
NDC effect is confirmed herewith also for systems with 
the energetic disorder. However, at extremely large 
energetic disorder (parameter a), the peak in the field 
dependence of the drift velocity disappears completely. 
The effect of energetic disorder becomes smaller at 
larger fields, as expected from the fact that in the limit 
F — ► oo the hopping rates do not depend on site energies. 



FIG. 12: Drift velocity v versus electric field F for Gaussian 
disorder in site energies. Disorder is characterized by a stan- 
dard deviation a of site energies from the reference energy. 
The localization length a — 0.15 R. 

by both Eq. Q and Eq. (|24|) . in accordance with the re- 
sult that the simulated NDC effect is weaker than the 
predicted one. We would like to emphasize that the nu- 
merically obtained NDC has exactly the origin predicted 
by Nguyen and Shklovskii, 15 consisting in spreading of 
the optimal trap into the region x < at finite fields 
and consequently in the increase of the trap volume with 
decreasing the field strength. The possibility of traps in 
the form of clusters instead of single chains of sites has 
been considered in Ref. We interpret our numerical 
result as a confirmation of that idea. 

A further decrease of the electric field F results in the 
washing out the empty region in the density of sites p(r) , 
as shown in Fig. UTTc) for F — F m . Finally, at small F 
the trap almost disappeared (Fig. [TTT d)). which points 
to a negligible role of the trapping effect in the regime of 
Ohmic conduction. 

Materials studied experimentally usually have disorder 
not only in the spatial distribution of localized states, 
but also in the site energies i 3 ' 5 ' 36 It is therefore necessary 
to check how stable the NDC effect is with respect to 
energetic disorder. In order to study the role of the 
energetic disorder for the NDC effect, we repeated the 
simulation in a system with a Gaussian distribution of 
site energies, characterized by the standard deviation a. 
Fig. [T2] shows that introducing a random energy for each 



VI. CONCLUSIONS 

Numerical studies of the field-dependent drift velocity 
of charge carriers in the hopping regime at high electric 
fields confirm the prediction of the existing analytical 
theorie s 14 ' 15 that the negative differential conductivity is 
inherent for this transport mode. However, the shape 
of the field dependence on the current density obtained 
numerically differs essentially from the one predicted so 
far^ The analytical theory has been improved to give a 
much better agreement with the numerical results. In the 
limit of the infinitely high electric fields, the predictions 
of the analytical theory of Nguyen and Shklovskii are 
to much extent confirmed by our straightforward Monte 
Carlo simulations. 
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